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To recover the flow information encoded in travel-time data of time-distance helioseismol- 
ogy, accurate forward modeling and a robust inversion of the travel times are required. We 
accomplish this using three-dimensional finite-frequency travel-time sensitivity kernels for 
flows along with a 2+1 dimensional (2+ ID) optimally localized averaging (OLA) inversion 
scheme. Travel times are measured by ridge filtering MDI full-disk Doppler data and the 
corresponding Born sensitivity kernels are computed for these particular travel times. We 
also utilize the full noise covariance properties of the travel times which allow us to accurately 
estimate the errors for all inversions. The whole procedure is thus fully consistent. Due to 
ridge filtering, the kernel functions separate in the horizontal and vertical directions, motivat- 
ing our choice of a 2+ ID inversion implementation. The inversion procedure also minimizes 
cross-talk effects among the three flow components, and the averaging kernels resulting from 
the inversion show very small amounts of cross-talk. We obtain three-dimensional maps of 
vector solar flows in the quiet Sun at spatial resolutions of 7 — 10 Mm using generally 24 h 
of data. For all of the flow maps we provide averaging kernels and the noise estimates. We 
present examples to test the inferred flows, such as a comparison with Doppler data, in 
which we find a correlation of 0.9. We also present results for quiet- Sun supergranular flows 
at different depths in the upper convection zone. Our estimation of the vertical velocity 
shows good qualitative agreement with the horizontal vector flows. We also show vertical 
flows measured solely from f-mode travel times. In addition, we demonstrate how to directly 
invert for the horizontal divergence and flow vorticity. We finally study inferred flow-map 
correlations at different depths and find a rapid decrease in this correlation with depth, 
consistent with other recent local helioseismic analyses. 

I. INTRODUCTION 

Time-distance helioseismology [1] is a set of tools that measures and interprets the travel times 
of seismic waves propagating from one point on the solar surface to any other point. It has been 
shown that these travel times contain information about solar flows [21 El [H El El among others] . 
This paper focuses on the inversion of travel times to obtain high spatial resolution maps of near- 
surface vector flows in quiet-Sun regions. What is unique to this study is that it is the first 
fully consistent inversion in time-distance helioseismology. The consistency is described by several 
factors: (1) we measure travel times with the same definition with which the travel-time sensitivity 
kernels are computed; (2) we use three-dimensional (3D) finite- frequency Born sensitivity kernels 
which are necessary to detect flow structures that have spatial scales on the order of the mode 
wavelength - this is the regime where the commonly used ray approximation fails for example] ; 

(3) we use the full noise covariance properties of the travel times as an ingredient in the inversion; 

(4) the inversion procedure we choose to implement is 'optimal', in that it simultaneously achieves 
the best possible spatial resolution while minimizing the magnification of the errors. Furthermore, 
the regularization is carried out in both the horizontal and vertical directions. 

We have developed a novel two plus one dimensional (2+ ID) inversion scheme based on the well- 
known subtractive optimally localized averages (SOLA) technique [HE]. The inversion procedure 
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explicitly minimizes the cross-talk effects among the three flow components by imposing constraints 
on the averaging kernels. An important aspect that we introduce to this procedure is that of ridge- 
filtered travel-time measurements, whereby only wave packets of a particular radial order are used. 
Sets of point-to-annulus travel times are then computed for the f, pi, P2, P3, and P4 ridges. This 
is quite different than the usual phase-speed filtering implemented for travel-time measurements. 
Subsequently, the sensitivity kernels are also computed as corresponding ridge-filtered quantities 
to match the travel times. The 2+ ID inversion is motivated by the observation that to a very 
good approximation the 3D Born ridge-filtered sensitivity kernels separate into the product of a 
2D horizontal function and a ID function in depth. 

Together, these ingredients allow us to infer all three components of the vector flow in the 
near-surface layers of the quiet-Sun convection zone. In previous work, we determined for the 
first time the maximum amplitude of flows which can be reliably recovered from a linear model of 
the travel-time perturbations [TO] . The supergranular and other quiet-Sun flows (with velocities 
< 400 ms- 1 ) that we detect in this study fall within this range, giving us confidence in the reliability 
of the method. The horizontal spatial resolution of the inversion presented here (^7—10 Mm, 
depending on the observation time, depth, etc.) is on the order of, or even in some cases below, 
the wavelength of the waves used in the analysis (typically 5 — 20 Mm, depending on the dominant 
modes). In addition, we perform a direct measurement of the vertical component of the flow 
(without simply invoking mass conservation) and with confidence that the cross talk between the 
horizontal and vertical components has been minimized that minimizes. Interestingly, we find that 
we can determine the vertical velocity even from f-mode travel times. We also directly invert for 
the horizontal divergence of the flow as well as the vertical component of the flow vorticity. 

As the main aim of this paper is to develop the inversion procedure and to perform tests of it 
on real solar data, we postpone the main interpretation of the results to a future publication. We 
typically show quiet-Sun flows that have been obtained from 24 h of data, in order to maximize 
the signal-to-noise ratio from the supergranulation signal [IT] . In all of the inferred flow maps we 
provide estimates of the noise and the spatial resolution. 

The paper is organized as follows: In Section [IT] we describe the data and the ridge-filtered 
travel-time measurements. That is followed by a brief discussion of the forward modeling, i.e., the 
computation of sensitivity kernels and the noise covariances that are consistent with the travel 
times. Since the multi-step inversion procedure is somewhat complicated, we provide an overview 
in Section |IV[ followed by two sections that discuss in detail the 2D and ID parts with plenty 
of example calculations. Three-dimensional averaging kernels at different depths are presented in 
Section [VII, and finally the various flow maps are presented, tested, and discussed in Section VIII| 
We end with a summary of the results and a discussion of current and future work. 



II. DATA AND RIDGE-FILTERED TRAVEL-TIME MEASUREMENTS 

For this study we use dopplergram data from the Michelson Doppler Imager [12] on board 
SOHO which are full-disk images of 0.12° spatial sampling (~ 1.46 Mm) and 1 minute cadence. 
The region of interest is an area centered on NOAA AR 9787 observed between 20-28 January 
2002 and tracked and remapped courtesy of T. L. Duvall Jr. [28] The size of the full set of velocity 
data cubes is 512 x 512 x 1440 x 9 (two spatial dimensions, 1440 minutes, 9 days). We only use 
the middle seven days for the results shown here. These data are ideal for helioseismic analysis as 
there is a sunspot that is large, isolated, and quite stable as it traverses the disk. For our purposes, 
there are also large regions of relatively quiet Sun in these maps where we will focus our analysis. 

We denote a Doppler velocity cube as 0(r,t), where 



r (x,y) 



(1) 
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is the horizontal coordinate, x and y are the east-west and north-south directions, respectively, 
t is time, and we work in a Cartesian geometry. We filter the data by multiplying the Fourier 
transform of the data cube by a filter F n {k,u) which selects all modes with the same radial order 
n, and removes all others. We call this ridge filtering. The ridge-filtered data <E> n are then given by 

^ n {k,u) = F n {k,u)<l>{k,u\ (2) 

where k is the horizontal wavevector and u is the angular frequency. The ridges we filter and retain 
for this study are the surface-gravity wave (f-mode) ridge and the first four acoustic (p-mode) ridges. 
Throughout the text we carry the index n which takes the possible values n — {f , Pi, P2, P3, P4}- 
We refer to these as mode ridges. It is important to note that this type of filtering is different from 
the phase-speed filtering that is typically done prior to any time-distance analysis. 

The cross-covariance functions are computed from <& n (fc,u;) for three different point-to-annulus 
geometries, denoted by 'oi' (out minus in), 'we' (west minus east), and 'ns' (north minus south). 
The 'oi' covariances are measured using the wave signal at a given point and the wave signal 
averaged over a concentric annulus of radius A. The 'we' ('ns') quantities use the central wave 
signal along with the wave signal averaged over the annulus but weighted by cos 9 (sin 9) , where 9 
is the angle between the x direction and each point on the annulus. This particular procedure was 
introduced in [9] and is similar to what is usually done in standard time-distance measurements [3] . 
The temporal cross-covariance functions are computed for 20 different annulus radii (A = 1.46 Mm 
to 29.2 Mm, incremented by dA = 1.46 Mm) for each mode ridge and measurement type, which 
yields a set of functions C a,n (r, t; A), where a = {oi, we, ns}. Flows introduce asymmetries in time 
lag t in the cross-covariance functions. 

Travel times are then measured according to the procedure developed by [HI [13]. This method 
establishes a linear relationship between the travel-time perturbations r and the cross covariances, 
given by 



W n (A, t) C Q ' n (r, t; A) - C ref ' n (r, t; A) 



(3) 



where the sum over t denotes a discreet sum over the observation time T (for this work T — 1 day 
in most cases), h t = 1 min. is the temporal sampling rate of the data, W are weight functions, 
and C ref is a reference cross covariance symmetric in time lag, which we choose to be the cross 
covariance computed from our model power spectrum. Thi s model is tuned to match the observed 
power spectrum of each individual ridge (see Section III). C ref is the inverse Fourier transform 
of the model power spectrum. Full details about the W function and equation ^ can be found 
in [llj . The resulting three different types of travel-time measurements are constructed to have 
sensitivity to different flow geometries. The 'oi' travel times are highly sensitive to horizontal flow 
divergence, whereas the 'we' and 'ns' travel times give information about directional flows. 



III. FORWARD MODELING 

We now briefly discuss the steps we have taken to model the ridge-filtered travel-time measure- 
ments described in the previous section, as well as the measurement noise properties. Full details 
of forward modeling in time-distance helioseismology can be found in [13], [TT] , and [T3j . 
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FIG. 1: Examples of various sensitivity kernels K a ' f (r, z; A) used in this study and defined in equation Q. 
They have been obtained by computing weighted azimuthal averages of the point-to-point kernels for the a 
point-to-annulus geometries. The kernels in the left column give the sensitivity of 'oi' f-mode travel times 
to u x , and the kernels in the middle (right) column give the sensitivity of the 'ns' ('we') f-mode travel times 
to u y (u z ). The top row are 2D kernels after integrating the 3D kernels over depth. The bottom row shows 
depth slices along y = of the 3D kernels. The type a and mode ridge n of the kernels in each column is 
indicated in the bottom panels. For each case, A = 10.2 Mm. The gray scale in the bottom row has been 
truncated to 75% of kernel maximum. 



A. Travel-time sensitivity kernels for ridge filtering 

We consider travel-time measurements of type a — {oi, we, ns} for each distance A and for each 
mode ridge n as described in the previous section. The travel-time perturbations are related to 
the small-amplitude flows through a linear relation: 

r n Q (r,;A) = h^.h z K a ' n (rj — r*j, z; A) • u(rj, z) + A/^' n (r*j; A), (4) 

3,z 

where K denotes the three-dimensional travel-time sensitivity kernel, u is the real vector flow in 
the Sun, Af tt represents the noise in the travel times (tt), and the sum over the vertical coordinate 
z denotes a discreet sum (note z = at the surface and is negative inside the Sun). The kernels 
are computed so that the horizontal grid spacing, h x = h y = 1.46 Mm where h% = h x h y , matches 
that of the travel-time measurements. The vertical grid (with spacing h z ) is taken from the model 
on which the kernels are computed. The kernel K is computed in the first Born approximation 
[T3] . We start with the point-to-point kernels for flows derived in [Mj, computed for the f, pi, p2, 
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P3, and p4 ridges, whose input power spectra have been tuned to match the ridge- filtered observed 
power spectra. For each ridge, 20 kernels are computed, one for each A. The 3D point-to-point 
kernels are then azimuthalfy averaged according to the three different point-to-annulus weighting 
geometries a used for the travel-time measurements. The total number of point-to-annulus kernels 
thus obtained is 300, each kernel having 3 components. A few examples of f-mode kernels after 
weighted azimuthal averaging are shown in Figure [I] for V{ — 0. From left to right, the columns 
show kernels that give the sensitivity to u x , u yi and u Zl respectively. The top row shows the 
resulting 2D kernel after integration over depth. 

As mentioned in the introduction, the motivation behind the inversion we choose to perform is 
that the ridge-filtered kernels to a very good approximation separate into a horizontal 2D function 
of r times a ID function of depth z. In other words, dropping the labels a, n, and A for the 
moment, the ith component of the sensitivity kernel defined in equation Q may be written as the 
product 

Ki(r,z)^fi(r) gi (z), (5) 

where 



h r h z J2j,z K i( r 3i z ) ^h 2 r h z J2j, z Ki(rj 

In Figure [2] we demonstrate the separability for an example p2 and 'we' averaged kernel. Figure [2^i 
shows a slice of the K^ e ' P2 kernel at about 100 km above the photospere, and in (b) we show the 
same slice, but of the function obtained by equation In the bottom two panels depth slices 
of these two kernels are shown for comparison. The match between the kernels is quite good. 
We have also studied the point-to-point kernels in this way. In general, the weighted azimuthally 
averaged kernels (the ones actually used in the inversion) separate 'better' than the point-to-point 
ones, because much of the small-scale structure is averaged away. All of the other kernels for ridge 
filtering that we have studied separate this way to a very good approximation, giving us confidence 
in the type of inversion we choose to employ. 



B. Noise covariance matrix 

Travel times contain a significant amount of realization noise and a good understanding of these 
noise properties allows us to assign accurate errors to the flow estimates. It has been shown in 
previous time-distance inversions that it is important to take into account the noise covariance 
matrix [151 EE E] • [H] showed in detail how to compute model noise covariances of travel times. 

We assume the solar oscillations are stationary and spatially homogeneous since we are re- 
stricting our study to the quiet Sun. We further assume that the noise between different ridge 
measurements n and n' is uncorrelated. This approximation is acceptable for the type of ridge 
filtering that is used in this study. The covariance matrix A of the noise components Af tt from 
equation Q is given by 

Af(n - rj; A, A') = Cov [tf£ n (n; A), A/^(r j; A')] . (7) 

This quantity has units of s 2 , and is computed according to equation (28) in [llj. Example plots 
for the case when n = f and for A = 5 Mm were shown in [9J using the same model. Similar 
features are seen for all mode ridges and distances considered here. Note that the covariance 
matrix elements of A scale with the observation time T as T _1 . 
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FIG. 2: Separability of a sensitivity kernel for p 2 and for ridge filtering. This particular kernel, K^ e ' P2 , gives 
the sensitivity of 'we' p2-mode travel times to u x for A = 23.4 Mm . (a) Horizontal cut through the kernel 
at a height of about 100 km above the photosphere, denoted by the dashed line in the panel below, (b) 
Horizontal cut through the quantity obtained as the product of the two (normalized) functions, f(r) and 
g(z), computed by integrating the kernel in (a) horizontally and over depth, according to equations ( |5|6| ). 
Panel (c) is a depth slice along y = of the original kernel from (a), and panel (d) is a depth slices along 
y = of the kernel in (b). The match between the left and right columns is quite good. The gray scale has 
been truncated to 40% of the maximum value for ease of comparison. 



IV. BASIC STRATEGY OF THE 2+1 DIMENSIONAL SUBTRACTIVE OPTIMALLY 

LOCALIZED AVERAGES INVERSION 

The problem we wish to solve is to estimate, for example, u x (r,z) in equation Q, given the 
travel-time measurements, the sensitivity kernels, and the noise covariance matrix. We carry this 
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out using a 2+ ID SOLA inversion procedure. We formulate the problem in terms of the component 
u x only for notational simplicity, noting that the procedure for finding all other flow components 
is completely equivalent. We first discuss the basic idea of this method, and then in the following 
two sections separately describe the 2D and ID parts in more detail. 

An OLA-type inversion for our purposes seeks a way to combine the sensitivity kernels to find 
a three-dimensional averaging kernel that is, roughly, the shape of a ball, centered horizontally 
about the origin r = and vertically about the target depth z t located somewhere in the solar 
interior. The averaging kernel will be found from the 2+ ID inversion and is defined by 

v x (r; z t ) = h 2 r h z 2 + 1D /C^ (rj - r, z\ z t ) • u(rj,z) + noise , (8) 

3,z 

where v x is an estimate of the real solar flow u x (as is the case in the rest of the paper) and the 
noise term will be specifically quantified below. To clarify the notation, sensitivity kernels are 
written as K and averaging kernels are written as JC. Any superscript to the left of the averaging 
kernel denotes by which type of inversion it was computed, for example, 1D /C, 2D /C, and 2+1D /C. The 
superscripts to the right indicate for which component of u and ridge n the averaging kernel is 
computed. 

It is important to observe from equation Q that the y and z components of the averaging kernel 
should be zero so that the flow estimate v x is not contaminated by any cross talk from u y and u z . 
As seen in the next section, the 2D inversion attempts to accomplish this by constraining the spatial 
integral of K y and K z to be zero. Ideally, one would want the x component of the averaging kernel 
to be a delta function; however, noise and a finite set of travel times inhibit this. Nonetheless, if 
an acceptable averaging kernel is found, then the travel times can be properly averaged to give an 
estimate of the local flow v x « u x in which we are interested. 

Since the problem essentially separates into a 2D and a ID problem because of ridge filtering, 
the three-dimensional averaging kernel 2 + 1D /C^ from equation ([§]) is derived in two main steps. 
The first step is to compute the 2D (horizontal) component of the averaging kernel, 2D /C(r), such 
that its x components is highly peaked about the point r = 0, and the other components are 
zero. This typically involves trying to match the x component to a target function that is a 2D 
Gaussian function in horizontal coordinates r = (x,y). The inversion coefficients, or weights w. 
that accomplish this averaging of the sensitivity kernels, are then used to combine the travel times 
in such a way that the estimated flow for any ridge measurements n is 

<(r; A) = £ ™ a ' n (rj ~ r; A)r£(r i; A). (9) 

The resulting flow map v™ is an average of the real flow u x over the depth that the dominant 
modes of ridge n probe. An intermediate step is to combine these maps over all distances A by an 
averaging procedure based on the correlated noise in the measurements. A set of weightings 7 is 
computed, described in detail in Section [VAj such that the distance- averaged flows are given by 

<(r) = 5>?<(r;A,-), (10) 



where 7j is defined in equation (21) and the sum j runs over all A used in the problem. 

The second main step is to obtain localization of the 3D averaging kernel in the vertical direction 
about zt by combining separate 2D ridge measurements. A ID inversion in depth is thus performed 
which seeks a new set of inversion coefficients c n . These new coefficients combine the 2D flow maps 



in equation (10) in such a way that the final estimate of the flow around zt is given by 

v x {r\ z t ) = V c n (zt)v%(r) « u x (r; z t ), (11) 
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FIG. 3: Example distance weightings from a 2D f-mode inversion for u x using a target resolution of 11 Mm, 



as described in Section V A The observation time is 24 h. The open squares are the weights at each distance 
obtained by assuming that the the estimated flows v f x (r; A) are uncorrelated for different A (eq. 17 ). The 



points with filled circles are the weights taking the correlations properly into account (eq. [21]). The noise 
(7 for each case indicated in the legend. 



where u x {r\ z t ) denotes the real flow at a particular depth z t . The whole procedure is in principle 
carried out to estimate each flow component (u x , u y ,u z ) at many different target depths z t to infer 
the vector flow u{r,z) throughout a desired interior region. We now describe in more detail how 
the 2D and ID inversion weights are computed in the following two sections. 

V. 2D HORIZONTAL INVERSION 



Based on the separability of the sensitivity kernels discussed in Section III, the 2D inversion is 
formulated to solve equation Q for v x (r) using the 2D depth-integrated kernel 



K a ' n (r; A) = h z J2 K a ' n (r, z\ A), 



(12) 



which we compute for all a, n, and A available. 

We can define the averaging kernel that we wish to find from the 2D inversion, 2B X, Ux,n , by 
plugging equation Q into equation Q to obtain 

<(r; A) = h%J2 2D K u *' n (ri - r, z; A) • u(r u z) + £ w a ' n {r j - r; A)M% n (r r , A), (13) 



where 



3 /C^' n (r,z;A) = J^«j a ' n (r j ;A)K a ' n (r-r j ,«;A). 



(14) 
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FIG. 4: Two-dimensional (top panels) and three-dimensional (bottom panels) averaging kernel for a 2D 
inversion for u x using only f modes. The panels in the top row show the integral over depth of the corre- 
sponding 3D kernel below it. The bottom row shows depth slices of the three components of the 3D kernel 
along the y = line. The red line outlines the FWHM= 7.3 Mm of the 2D Gaussian target function used 
in the inversion (see equation [l5]). The overplotted color contours, which are also marked on the colorbar 
for reference, denote the following: (white) is the half maximum of the kernel, (black) is the negative of the 
half maximum, (blue) and (green) denote ±5% of the maximum value of the kernel, respectively. Estimates 
of the noise from this and all other inversions are given in Table [TJ 



This shows explicitly that the averaging kernel gives only an estimate v™ of the flow that is some 
average of the real flows over some depth. 

The details for obtaining the inversion weights w for a 2D SOLA inversion for flows were 
presented in [9j, however, in that work only one mode ridge, one distance A, and two components 
of the sensitivity kernels (K x and K y ) were utilized. The generalization to our present case is 
straightforward. To summarize this procedure, we first prescribe a target function 2r>/ T Ux that 
we wish the averaging kernel to resemble. It is a vector-valued function, chosen such that the x 
component is typically a 2D Gaussian in r with dispersion cr, and the other components are zero: 



-T-(r)=^^-^,0,0j, (15) 

where r = \\r\ \ is the 2D vector norm. The horizontal integral of the target function is normalized to 
one. The full-width at half-maximum (FWHM = 2aV2 In 2) of the target function is a measure of 
the resolution of the inversion if the averaging kernel matches it well. For the sake of completeness, 



in an inversion for the jth component of and denoting the gaussian function in equation (15) as 
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G(r), the ith component of the target function is 2D 7~^ = G(r)eiSij, where is the unit vector in 
the zth direction and 6 is the Kronecker delta function. 

Two quantities, one which measures the mismatch between the averaging kernel and target 
function, the other which measures the noise propagation, are computed. Let the noise in the 



inversion for u x for ridge n be denoted by J\f Ux,n ; it is specifically defined in Section |V A| A 
minimization (with respect to the inversion weights) is carried out according to 



mm 

w 



J2\\ 2n K Ux > n (ri) - 2B T Ux (n)\\ 2 + {N u ^ n f , (16) 



where (3 is some regularization parameter that we choose typically to be quite small [9J. A large 
matrix is then regularized and inverted for each value of the trade-off parameter, which results in a 
unique set of weights at each point in this parameter space. We choose an 'optimal' set of weights 
w from examining the trade-off curve (L curve) as discussed in [9J , such that the averaging kernel 
matches closely the target function. The 2D averaging kernel is constructed by convolution of the 



sensitivity kernels with the weights according to equation (14). 

Since the inversion in this example is carried out for u™, it is important that 2D /Cy x ' n and 2r> fCz x ' n 
be as close to zero as possible, to minimize the 'cross talk' among all of the components. This is 
achieved in practice by constraining the total spatial integrals of the y and z components to be 
zero, although in practice there is usually some structure present even with this constraint. We 
are in the process of exploring other effective constraints. When a well-localized averaging kernel 
is found for each ridge n, the weights are then suitable to be used to average the travel times to 
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give an estimate of the flow v™(r) using equation 



A. Combining all of the distances A 

Throughout this inversion procedure, it is necessary to combine the quantities we obtain for 
different annulus radii A, such as the estimated flow maps v™{r\ A). This is done by weighting each 
distance by appropriate weights. One simple way of achieving this, typically used in helioseismology, 
is to assume that the noise in each measurement is independent and uncorrelated. Then the 
standard deviation in the estimated flow maps is used to determine the contribution of the errors 
at each distance. We denote the standard deviation of a set of flows v™ for each distance Aj and 
mode ridge n as cr™ c > n . The 'unc' superscript emphasizes the assumption of uncorrelated data. 
Finding the minimum variance of this set then gives a weighting factor, 77™, according to 

(lK nc ' n ) 

V?= V / J N2 , (17) 

E20 / -1 / unc.Tl \ 
3=1 \}l a i ) 

where the sum in the denominator runs over all 20 distances used in this problem. We carry out a 
2D inversion as described above for u x using f modes in a region of quiet Sun (the same region used 
m Section [Villi ). The weights obtained from the estimated flows using equation (h 7b are plotted in 



Figure [3] as open squares. For very small distances where the noise is very high, the weights are 
zero. From the total variance we obtain a noise estimation given by 



a unc,f 



-1/2 



(18) 



which for this particular example is 14 ms _1 . 

However, we know that the values of the flows at different A are correlated quite strongly due 
to noise |11| , and so we choose to average them in a way that takes these correlations into account. 
A covariance matrix C n of the noise in the individual flow measurements at distances A and A / of 
ridge n is computed using the 2D inversion weights as 

C n (A, A) = ™ a ' n (rf, A)Af (r* - r j5 A, A V ,n fo; A), (19) 

i,j,a,f3 

where A is the covariance matrix of the noise in the travel times (equation [7]), and w are the 2D 
inversion weights. Note that the matrix C has units of length 2 time -2 . The final measurement of 
any general quantity q for each mode ridge is then obtained by averaging the 20 distances we use 
according to (for example, see [IE]): 

20 

q n (r) = Y,^q n (r;A j ), (20) 

3 = 1 

where the set of weightings 7" is given by 



„ [Yf =1 C-\A u A 3 )) 

it = 9 — (21) 
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FIG. 6: Two-dimensional (top panels) and depth slices of three-dimensional (bottom panels) averaging 
kernel for a 2D inversion for u z using only f modes. The colors are as for Figure [4j 



and the variance (cr 00 "' 71 ) 2 for the correlated case is 

(o-».») 2 = [^(C-^Ai.A,-)) 



(22) 



We can now identify JV Ux,n = cr corr ' n , the noise in a measurement of v™, introduced in equation (161). 
A set of weights for the f-mode case 7? obtained this way are plotted as the filled circles in Figure 3(to 
compare with the uncorrelated case. What is interesting to note is that for some distances a negative 
contribution is needed to average the data properly, which is never the case for the uncorrelated 
data. The estimated noise is also lower than in the uncorrelated case. Equation (20) is quite 
general, and has been used to obtain the distance-averaged flow maps presented in Section VIII , 
as well as all of the averaging kernels shown in this paper. Quantities written without the distance 
argument A have been averaged this way. 



B. Averaging kernels from the 2D inversion 

Example averaging kernels after combining all distances for a 2D inversion are shown in Fig- 
ures |4][7| In all figures, the top row shows the 2D averaging kernel obtained by integrating over 
depth the adjacent 3D averaging kernel in the bottom row. The bottom panels show depth slices 
along y = 0. The inherent noise from the inversion corresponding to each figure is given in Table [IJ 
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FIG. 7: Two-dimensional (top panels) and depth slices of three-dimensional (bottom panels) avera^ 
kernel for a 2D inversion for u z using only pi modes. The colors are as for Figure [4j 



;mg 



Figures [4] and [5] show averaging kernels for a 2D inversion for u x for the f and pi mode ridges, 
respectively, while Figures [6] and [7] are for an inversion for u z . The absence of any dominant cross 
talk is evident in the top panels of all figures, which are the 2D averaging kernels that come straight 
out of the 2D inversion. There is a completely negligible y component in all cases for the inversions 
for u x . The cross talk is slightly more pronounced for the x and y components of the kernels in 
the inversions for u Z: but is confined to the very near-surface region, typically above the depths 
which are significant for our inversion results. In general, these averaging kernels are quite good; 
the only structure from the 'off diagonal' terms are of the order of about 5% of the diagonal terms. 
It is important to study averaging kernels such as these to have an idea of what the inversion is 
actually accomplishing. Similar plots have been examined for all of the other mode ridges available 
and they exhibit similar features. 



C. Minimum variance 



The averaging kernels from the 2D inversion in Figures [4][7] are computed for each separate mode 
ridge n. One possible way to combine them over all ridges is to use a simple minimum variance 
treatment of the noise in the estimated flow component v™, as described in Section VA for the 
case of distance averaging. Because of ridge filtering, we consider the noise between mode ridges 
n and n f to be uncorrelated. Thus, any quantity q n that depends on the set of mode ridges can be 
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FIG. 8: Three-dimensional minimum variance averaging kernel for a 2D inversion for u x . The top panels 
are the 2D averaging kernels after integrating the 3D kernel over depth, and the bottom panels are slices at 
y = through the 3D kernel. The method to obtain this kernel is discussed in Section V C The resolution 
of the inversion is 7.3 Mm. The colors are as for Figure [4j 



averaged according to 



where the weights d n are given by 



n 



d n q u 



(23) 



cf 1 



(l/o- 



En (V^ 



(24) 



and a corr ' n , the correlated noise estimated from the v™ measurements, is denned in equation (22). 



In Figure [8] we show the three components of an averaging kernel K, Ux obtained by combining five 
kernels JC Ux,n from the minimum variance in v™ using the weights given in equation ( [24] ) . The 
minimum- variance weights for each ridge and noise for this figure are given in Table [IJ The noise 
level is low and the kernel is well localized horizontally, however, we clearly have no control over 
the depth at which one wishes to have sensitivity. 

In conclusion, the 2D inversion produces averaging functions (such as those shown in Figures [4j{8]) 
that are 'optimally' localized in the horizontal direction, but not in depth. This depth localization 
is accomplished by a ID inversion, to which we now turn. 
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the legend) for an inversion for u x . These are obtained according to equation (26). 
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FIG. 10: Example trade-off curves computed in the ID inversion for two target depths. The individual curves 
in each set of depths are obtained by varying the horizontal target function width. For the z t = — 1 Mm 
set, curves a, b, c, and d correspond to 5.8, 7.3, 8.7, and 10.2 Mm, respectively. For the z± — —3.5 Mm, a, 
b, c, and d corresond to 7.3, 8.7, 10.2, and 11.6 Mm, respectively. In this scheme, sets of inversion weights 
would be chosen at the points given by the circles. 



VI. ID SOLA DEPTH INVERSION 

Up to this point, the 2D inversion has provided averaging kernels which average the real solar 



flows to give estimates of the flows for each particular ridge, according to equation (13). If we 
assume that the real flows vary slowly in r over the horizontal extent of the averaging kernel, we 
can perform the summation in equation (Jl3]) over T{ to obtain 



<(r) *h z J2 1D K^(z)u x (z) + £ w^(r 3 - v)^^), (25) 



3,OL 
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where 

^K^iz) := h 2 r J2 2U K U x x ' n (n,z), (26) 



is the ID sensitivity kernel for an inversion for u x and mode ridge n. Recall that the horizontal 
integrals of 2D /C y and 2r> K, z are zero due to the constraint imposed in the 2D inversion; therefore, 



only the x component of the quantities remains in the right hand sides of equation (25) and (26). 

The five available ID sensitivity kernels are shown in Figure [9| The ID SOLA inversion seeks 
inversion coefficients c n (zt) that combine each ridge measurement about target depth zt, so that 



the final estimate of the flow using equation (25) is 
v x (r;z t ) = J2c n (zt)v n x (r) 

n 

= h z J2 10 ^ x (z; z t )u x (r, z)+Y^ c n {z t )w^(r 3 - r)N% n (r 3 \ 

z j,Oi^n 

where 

^(z;z t ) = J2c n (zt) 1D K x ' n (z) (27) 

n 

is the one-dimensional averaging kernel peaked about Zf. The ID coefficients c n are obtained in 
an analogous way to the 2D case. A target function is chosen which is typically a ID Gaussian in 
depth, centered about zt\ 

„-(z-z t ) 2 /2a 2 

™T x ^(z;z t ) = 7 =— . (28) 

A misfit quantity is constructed which measures the mismatch between the target and averaging 
functions: 

misfit 2 = h z [ 1D ^ (z; zt) - 1D T^ (z; z t )f . (29) 

Z 

In addition, a quantity which quantifies the error (noise) is included: 



error 



,2 \ A / n \ru x ,n\2 



(c n Af Ux > n Y . (30) 



A regularization parameter \i is introduced, and a minimization procedure with respect to the 
weights c n is carried out as 

min [misfit 2 + terror 2 ] . (31) 

Computing the minimization results in a system of linear equations, which is solved by inverting 
a small matrix to obtain the coefficients for each value of the regularization parameter. Finally, 
we choose weights roughly in the 'elbow' of the trade-off curve upon visual inspection. Several 



example trade-off curves from this procedure are shown in Figure [10] for two target depths z t and 
different target widths. 



In Figure [TT] we provide examples of ID averaging kernels for an inversion for u x and for different 
target depths. Due to a limited mode set, there is a limited number of depths that can be targeted 
properly. There is also an obvious limit on the maximum depth with which we can probe with 
these modes. In addition, as with other helioseismology inversions, a surface component is present 
(see, e.g., [H]). 
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FIG. 11: ID averaging kernels 1B JC Ux (z;z t ) for an inversion for u x for different target depths. These are 
obtained according to equation (27). 




FIG. 12: Vertical slice through a 3D averaging kernel 2+1D /Q a; after performing a 2+1D inversion for u x for 
a target depth of about —1 Mm below the solar surface. The left (middle) column is a slice along the y = 
(x — 0) line. The red contours show the half maximum value of the 3D target function, white contours show 
50% of the maximum value of the averaging kernel, black contours show —50% of the maximum value of 
the averaging kernel, and the blue (green) contours denote ±5% of the maximum. In the rightmost panel, 
the red line is the ID target function, and the black solid line is the ID averaging kernel. The ID weights 
c n and the noise estimate are given in Table [l| 
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FIG. 13: Vertical slices through an averaging kernel after performing a 2+1D inversion for u x for target 
depth zt ~ —3.5 Mm. The right panel shows the ID averaging kernel and target function. The colors are 



the same as for Figure 12 
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FIG. 14: Vertical slices through an averaging kernel after performing a 2+ ID inversion for u x for target 
depth zt ~ —7.5 Mm. The right panel shows the ID averaging kernel and target function. The colors are 



the same as for Figure 12 



VII. 3D AVERAGING KERNELS FROM THE 2+1D INVERSION 



We denote the final 3D averaging kernel produced from the 2+1D inversion for u x as 2 + ir> JC Ux . 
It has been defined in equation ([§]), and can be constructed from both sets of inversion coefficients 
(w, c) in terms of the original sensitivity kernels using equations (11) and (13): 



2+lDi 



K u *(r,z;z t ) := V c n (z t ) 2I) K u *> n (r, z 



n 



(32) 
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TABLE I: Values of the relevant quantities for some of the figures in this paper. Listed is the figure number, 
the type of inversion used to produce the figure, the ID inversion weights c n (or minimum variance weights 
gP), and the estimated noise associated with each measurement (1<j values are given). Note that the ID 
inversion weights do not apply for the 2D inversion. 
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:= ^^(zt^w^irfA^K^ir-r^z-Ai), (33) 

ijan 

where we emphasize that the weights w and c are obtained from a specific inversion for u x . Recall 
that the weights 7 are used to average the quantities over distance A. 

We now check to see if the final averaging kernels are as well localized as can be expected from 
the mode set used, which also justifies separating the problem into 2D and ID parts. Performing 
the full 2+ ID inversion for u x for different target depths produces example 3D averaging kernels 



such as those shown in Figures [12j [13j and [T4j Plotted in the left and center panels in each case 
are depth slices along the y = and x = lines of the x component of the kernel, 2+1D /Q x . The 
contour of the half-maximum value of the 3D target function is overplotted on the depth slices 
in red. The white contours show the half-maximum value of the 3D averaging kernel. Note that 
in principle, in a noiseless inversion with a large set of available modes, the white contours would 
match perfectly the red contours. The blue and green contour lines denote ±5% of the maximum 
value. Also provided in the rightmost panel in each figure for comparison are the ID target and 
averaging functions. For the shallowest target depth, instead of a simple ID gaussian we use a 
target function that goes to zero at z = 0. The ID inversion coefficients c n used to construct the 



kernels in Figures 12 - 14 are provided in Table [IJ The shallowest depth we can reach with these 



modes is about —1 Mm, the deepest about —8 Mm. Of course, the noise begins to increase quickly 
with depth. 



VIII. RESULTS WITH MDI DATA FOR QUIET-SUN FLOWS 

In the rest of the paper we provide example flow maps in the quiet Sun from the 2D and 2+1D 
inversion procedure discussed above. It is our main intention to demonstrate that the results 
obtained are sensible and consistent with what might be expected given the spatial resolution, 
observation time, and level of estimated noise. We note also that in order to study local flows, we 
typically remove a mean, large-scale, time- averaged flow from each retrieved map. 

The first simple test we perform is to compare the inferred flows from the 2D inversion to the 



20 



direct MDI Doppler data. To accomplish this, the three components of the inferred 2D vector flows 
are projected onto the line-of-sight vector at each pixel. We use 24 h of data from the seventh day 
(Jan. 26 2002) of the nine-day data set available. There is a sunspot in the center of the map 
which is centered vertically about the equator and at 30 degrees towards the western solar limb 
for this day. Figure [15] shows the comparison with the Doppler map from a 2D inversion using 
only f modes. The correlation between the inferred flows and the Doppler flows is 0.9 for pixels 
with less than |10| G of magnetic field. Also provided is a magnetogram to show the locations of a 
large sunspot, the surrounding plage, and the region of quiet Sun that is studied in all of the plots 
in the rest of the paper. The scatterplot of the two velocity maps shows that the magnetic field 
introduces an anomalous second component to the velocity field, reinforcing our reason to restrict 
all further analyses to the quiet Sun. That the inferred line-of-sight velocity is smaller than the 
Doppler velocity in not surprising. It is most likely attributable [5l [20] to the average depth over 
which the flows are measured. We cannot make a purely 'surface' measurement with the available 
mode set. 



Tests of the 2D inversion 



A further example of the 2D inversion is shown in Figure 16 Here we show horizontal flow maps 
inferred from inverting individual ridge travel times for f, pi, and p2 and for 24 h. The colorscale 
of these images is the horizontal divergence, obtained from a separate inversion as discussed below. 
In the f-mode map we see signatures of supergranulation with strong horizontal divergence in the 
centers of the supergranules and flow convergence at the cell boundaries. The outflow is generally 
in the 200 — 350 ms _1 range. The horizontal flows obtained from pi and p2 travel times are weaker, 
as is the supergranulation signature. Maps for the P3 and p4 ridges have also been studied, but 
the noise begins to increase quickly with these modes at this resolution. 

Obtaining the horizontal flow divergence from a direct inversion such as shown by the color 



scale in the plots in Figure 16 is conveniently done in an OLA inversion such as this one. All 
that is needed is to use a different 2D target function such that the quantity that is inverted for is 
simply the horizontal divergence. 

The function we use (which would replace the function defined in Eq. [15] in the 2D inversion) 



is 

/-7-divh 



where divh is a superscript label that denotes the horizontal divergence of the flow, Vh-^h- Another 
useful quantity we have studied is the vertical component of the flow vorticity, z • Vh x v^. We 
label the associated target function with the superscript vort^, which can be shown to be 

^)=(^e-^^ e -^,„Y ,35) 



An example map of the vertical vorticity obtained directly from a 2D inversion is shown in Figure [17 
We have studied plots of the divergence and vorticity at different resolutions and checked that the 
inverted quantities using these two types of targets and the direct numerical computation of these 
quantities using inferred horizontal vector flows agree reasonably well. The advantage of computing 
them directly from an inversion is that we are able to determine the noise and spatial resolution 
properly. 
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FIG. 15: Test of the 2D inversion. The top panels show a comparison of the inverted flows with the direct 
Doppler data for a map centered at 30 degrees west of disk center. The top left panel is the inferred vector 
flow map after projection onto the line-of-sight vector at each pixel for an inversion using f modes and 24 h 
of data. The resolution is about 7 Mm. The Doppler map is an average over 1 day of MDI full-disk data, 
smoothed with the 2D averaging kernel from the inversion and multiplied by the slope 0.61 of the best fit 
line through the scatterplot (panel below) to allow for direct comparison with the inverted map. The lower 
left panel is the 1-day averaged (truncated) MDI magnetogram with a sunspot visible in the middle. The 
units are in gauss. The white box outlines the quiet-Sun region analyzed for all of the plots in the rest of 
the paper. The white dashed line shows the location of the slice for the flows in Figure 24 The lower right 
panel shows a scatterplot of the velocity maps, where red (blue) dots represent pixels of less (more) than 
1 10 1 G. The correlation is 0.9 for the non-magnetic data. The dashed line is y = x, and the solid line a fit 
to the non-magnetic data (slope 0.61). The small vertical line represents the ±la (a = 17 ms" 1 ) noise in 
the flow estimation. 
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2D inversion using f modes. Noise in v x : 9 m/s. 
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FIG. 16: Horizontal flows and divergence from a 2D inversion using 24 hr of data for the f, pi, and p2 ridges. 
In each panel, the arrows denote the horizontal flows (obtained from inverting for v x and v y ) for 24 h of 
data, and the color scale is the horizontal flow divergence, obtained from a separate inversion as described 
in the text at the end of Section VIII A The x component of the 2D averaging kernel from the flow inversion 
is plotted in the box in the upper left of each panel and the FWHM, outlined by the circle, is 11.6 Mm. The 
x component of the 2D averaging kernel for the horizontal divergence is given by the quantity in the lower 
right box (see Eq. 34 ). Note the strong supergranular flows in the f-mode map which gradually weaken as 
the modes probe deeper layers of the convection zone. The correlation of the f-mode map with the pi-mode 
map is 0.88, and that between the f-mode and p2-mode maps is 0.35. The noise in v x is given for each panel. 
The noise in the horizontal divergence inversion is 10, 12, and 13 Ms -1 for the top, middle, and bottom 
panels, respectively. 
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2D inversion using f modes. Noise in v x : 17 m/s. Noise in v vortz : 5.9 Ms 1 
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FIG. 17: 2D inversions for the vertical vorticity and horizontal flows in the quiet Sun. The arrows denote 
the horizontal flows and the color scale is the vertical vorticity obtained from a separate inversion using 24 h 
or data. The x component of the averaging kernel from the inversion for the horizontal flows is given in the 
upper left box, and the x component of the averaging kernel for the vorticity inversion, which matches the 



target function from equation (35), is shown in the lower left box. 



B. Tests of the 2+ ID inversion 



A consistency test of the full 2+ ID inversion is shown in Figure [T8j We study flows in a quiet- 
Sun region obtained from two independent inversions using 24 h of data. Shown are the horizontal 
flows computed from inversions for v x and v y along with the line-of-sight magnetic field (colorscale). 
In the inversion corresponding to panel (a), we only invert pi travel times. We then attemp to 
target the pi averaging kernel by using all of the other available mode-ridge kernels except pi, i.e., 
f? P2 5 P3 5 and p4. Once a similar averaging kernel is found, we invert the corresponding f, p2, P3, 



and p4 travel times, and the resulting flows are shown in Figure [18p. The ID averaging kernels for 
each case are given in the panel on the bottom right. The maps are quite similar (correlation 0.82), 
and the differences could be due simply to the differences in the averaging kernels. In these maps 
the supergranule-scale flows are evident, and the magnetic field is concentrated at the boundaries of 
the supergranules as expected. A best fit through a scatterplot of the data taking into account the 
noise in both variables gives a slope of 0.83. The magnetic field does not introduce any anomalous 



component in the scatter as it did for the full-map study in Figure 15, confirming that this is a 
quiet-Sun region for our purposes. We have also used this test to see if we can recover an f-mode 
map by inverting the four available acoustic-mode travel-time sets. For as closely as we are able 
to match the averaging kernels, it is successful. The same conclusion can be drawn for the other 
possible cases when noise is not a limiting factor. 

In all of the plots studied so far, we have shown flows obtained from 1 day of travel times. In 



Figure 19, we compare horizontal flows at a depth of 1 Mm below the surface from inversions for 
different observation times. The panels show inversions for 6 h to 24 h in six hour intervals. Also 
shown is the horizontal divergence computed numerically (we have not yet computed inversions 
directly for horizontal divergence at depth). What is evident is that even with 6 h of data and a 
resolution of about 9 Mm, the noise level is reasonable and features are seen that have much in 
common with the 1 day map. The correlation between the 6 h and 24 h maps is still about 0.8. It 
is encouraging that the supergranulation signal at this depth is not dominated by noise for 6 h of 
data. 

We now compare horizontal flows at three different depths from the full 2+ ID inversion. Fig- 
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(b) 2+ ID inversion using f, D2, P3, and p4 modes. Noise in v x : 17 m/s 

1 L^, /// \ I • 1 1 1 1 1 ' 1 ' 1 ' 1 ' L 



- / / ^-^\y. / 



; . - 200 m/s 



100 



50 § 

bX) 



o cS 

.° 

a; 

-50 8) 

CO 



'-100 
1 100 



50 § 

bX) 



CP 




-200 -100 100 

Velocity using pi modes [m/s] (a) 



0.2 0.4 
Averaging kernel [Mm -1 ] 



FIG. 18: Test of the full 2+ ID inversion. Panel (a) shows the flows from a 2D inverison using only the 
Pi-mode ridge travel times. Panel (b) shows flows from a 2+ ID inversion using travel times from all other 
mode ridges (no pi modes) by attempting to target the same averaging kernel. The arrows denote the 
horizontal flows (obtained from inverting for v x and v y ) for 24 h of data, and the color scale is the truncated 
magnetic field from MDI. The x component of the 2D averaging kernel is plotted in the box in the upper 
left, and the FWHM= 11 Mm is outlined by the circle. The scatterplot shows the flows from panel (a) 
vs. (b), where the red (blue) dots are the values in the maps for which the magnetic field of the pixel is 
less (greater) than |5| G. The dashed line is the y = x line and the the solid line is the best fit through 
the non-magnetic data (slope= 0.83), taking account of the errors on each axis. The la error bars for each 
measurement are denoted by the cross in the lower right of the plot. The correlation in the scatter is 0.82. 
The bottom right panel shows the ID averaging kernel (solid line) for panel (a) and the ID averaging kernel 
(dashed line) for panel (b). The ID inversion weights for these two inversions are listed in Table [l[ 
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FIG. 19: Comparison of flows for different observation times T obtained from 2+1D inversions at a depth of 
1 Mm beneath the surface. The arrows are the inferred horizontal flows and the color scale is the horizontal 
divergence computed numerically. The total observation time increases from top to bottom in 6 h intervals, 
but the resolution remains fixed. The correlation of the 6 h map with the 24 h map is 0.78, while the 
correlation of the 12 h (18 h) map with the 1 day map is 0.9 (0.97). Note that the divergence scale is the 
same for each panel. 
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2+1D inversion at depth -1.0 Mm. Noise in v x : 10 m/s. 
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FIG. 20: Comparison of horizontal flows (arrows) at different depths from a 2+1D inversion. The background 
colorscale is the horizontal divergence to emphasize the various flow structures and obtained by numerical 
differentiation of the 2D flows. The x component of the 2D averaging kernels are shown in the upper left 
and all have FWHM= 11.6 Mm. The correlation of the x component of the flows at a depth of —1 Mm and 
—2.6 Mm is 0.77. The correlation of the x component of the flows at a depth of —1 Mm and —3.7 maps is 
0.33. The ID inversion coefficients for each depth inversion are provided in Table [l| 
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FIG. 21: Correlation with respect to the near-surface value of the inferred v x flows at different depths using 
24 h of data. The values represent the average over 5 days of the flows measured in the area of the quiet Sun 
used throughout this paper. The la error bars are plotted at each depth point, obtained from the standard 
deviation in the correlations over the 5 days. The 'near-surface value' for this case is at a depth of —1 Mm. 
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20 shows the flow field at depths of 1 Mm (top), 2.6 Mm (middle), and 3.7 Mm (bottom) below 



the surface using 24 h of data. The color scale is the horizontal divergence computed by numerical 
differentiation of v x and v y . The ID inversion weights for each map are given in Table [ij The flows 



at the different depths in Figure 20 are not too unlike the individual ridge flows shown in Figure 16 



and inspection of the ID inversion weights confirms that this should be the case. This figure also 
demonstrates that combining the maps with the ID inversion not only gives a good estimate of 
the target depth, but also generally lowers the noise levels. 

We have studied the correlation of maps of v x and v y such as those in Figure 20 at different 
depths with the near-surface map and averaged over 5 days. A plot of the results is provided in 



Figure 21, Each measurement is for 24 h of data, and the error bars are obtained by studying the 
variance in the correlation values. The correlation steadily decreases as we go deeper, and seems 
to disappear at about 5 Mm below the surface. However, the noise levels at these depths are quite 
large and we can draw no other specific conclusions at this time. This is consistent with recent 
studies on realistic numerical simulations using time-distance helioseismology [21 J and helioseismic 
holography [22]. In fact, the authors in [22] note that ". . . supergranule-sized flows are essentially 
undetectable using current methods below depths around 5 Mm ..." using 24 h of data or less. 
We confirm this conclusion here, and note that similar results have also been found with direct 
modeling techniques [23] . 

We have also studied the day-to-day correlation of the v x and v y maps at various depths. If 
we were predominantly measuring noise, there would be no significant correlation from one 24 h 
period to the next. Computing an average day-to-day correlation over seven days of data for the 
— 1 Mm depth maps gives a value of about 0.4. For the —3.7 Mm depth, we find a 0.26 correlation, 
and at a depth of —6 Mm, about a 0.1 correlation. This again demonstrates that there is plenty of 
near-surface flow signal when 24 h averages are studied, presumably due to supergranulation [TT] . 
which then quickly decreases with depth. 
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C. Vertical flows 



It has proven difficult in helioseismology to accurately measure the vertical component of the 
velocity near the surface due in part to its small magnitude compared to the other components. 
In fact, in many helioseismic inversions for flows, v z is approximately obtained by computing the 
horizontal flow component and then invoking mass conservation from the continuity equation (see 
[21] for an example in ring-diagram analysis). Another source of difficulty in these measurements 
has been associated with cross-talk effects, whereby the inversion (or sensitivity kernel) becomes 
insensitive to differences between upflows and convergence, and downflows and divergence [21]. 
These inversions, usually based on the ray approximation, have no obvious means of constraining 
the cross talk. Since we have available Born sensitivity kernels for v z , and an inversion procedure 
which measures each flow component while minimizing the cross talk with the others, we can 
obtain vertical flows directly and with the assurance that they are relatively independent from the 
horizontal measurements. This is clearly demonstrated in the averaging kernels of Figures [6] and 
[7J We note that we have so far only tested the 2D inversion for v z \ thus, the maps shown here are 
for individual ridge measurements. 

There tends to be much more relative noise in the measurements of vertical velocity, and there- 
fore in Figure 22 we show the vertical component of the velocity as the color scale averaged over 
2 days from a 2D inversion (the noise goes as T -1 / 2 , where T is the observation time). The top 
panel of Figure 22 is for the f-mode ridge, the middle panel for pi, and the bottom panel for 
P2- Also shown are the corresponding horizontal flows given by the arrows. One generally sees a 
good correspondence in all maps between the vertical upflows and horizontal outflows, as well as 
between downflows and horizontal inflows. Analysis of many similar maps show that the speeds of 
the vertical flows in the center of supergranules near the surface are on average about 15 — 20 % of 
the speeds of the horizontal outflow in the supergranules, slightly higher than recent observations 
might suggest [25] . 

To understand if the inferred vertical flows at these depths for p2 are reasonable, we compare 
them with maps of the horizontal divergence, obtained from a separate and independent inversion of 
the this quantity as explained in Section VIII A The vertical component of the flow and horizontal 
divergence are proportional if one writes down an approximate continuity equation whereby one 
neglects the horizontal variations in the density and the vertical gradient of the vertical flow. 
The scaling factor is the density scale height. In Figure 23 we show a scatterplot of v z 2 against 
the horizontal divergence inferred from inverting p2 travel times for the same region of the Sun 



as in Figure 22, The correlation coefficient is 0.62. The noise in v z is 9 ms _1 and 9 Ms -1 for 
the divergence measurement. The slope of the best fit line, using the noise information in both 
variables, gives a value of about 2.3 Mm. This value is in the range of the density scale height for 
the implied depth range of these vertical flows. We have also studied the correlation of vertical 
flows maps with horizontal divergence maps for the f-mode and pi-mode cases. The values are 
always in the range of 0.6-0.7. 

Another interesting question is how well the the near-surface vertical flows are correlated with 
deeper vertical flows. Since we have so far only implemented the 2D inversion scheme for vertical 
flows, we take different mode ridges as a proxy for depth. We correlate the 24 h f-mode v z map 
with the pi and p2 maps, average over seven days, and find correlations of about 0.3 and 0.2, 
respectively. In addition, as was described previously for the horizontal component, we have also 
studied the day-to-day correlations of the vertical flows averaged over seven days of data. The 
average day-to-day f-mode map correlation is 0.15, 0.2 for Vz 1 , and 0.15 for v z 2 . This demonstrates 
again that the v™ inversions are measuring long-lived flow structures and not just noise. 

Finally, in Figure 24 we show the culmination of our main results. It is a slice in depth through 
the horizontal divergence with overplotted v z information. The inversion to obtain these flows used 
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2D inversion using f modes. Noise in v x : 6.4 m/s. Noise in v z : 8.3 m/s. 
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FIG. 22: Vertical flows (color scale) and horizontal flows (arrows) in the quiet Sun for a 2D inversion using f 
modes (top), pi modes (middle), and P2 modes (bottom). These flows were obtained using 48 h of data. A 
positive vertical velocity means an upflow. The z component of the 2D averaging kernel from the inversion 
for u z is given by the quantity in the box in the lower right. The noise for all of the measurements is 
indicated. The correlation of these particular pi and P2 vertical flows with the f-mode map is about 0.3 and 



0.2, respectively. The region of the Sun used here is outlined by the white box in Figure [15] The color scale 
is the same in each panel to ease comparison. 
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FIG. 23: Vertical velocity versus horizontal divergence for flows using p2 modes measured over 2 days for 
the same region of the Sun as shown in Figure 22 The vertical flows and the horizontal divergence were 



computed from two separate inversions. The correlation coefficient is 0.62. The line shows a best fit through 
the scatter taking into account the noise on each axis, which is 9 ms _1 for v z and 9 Ms -1 for the divergence, 
indicated by the cross in the upper left of the figure. The slope of the best-fit line is a rough proxy for the 
density scale height at the implied depth. 



72 h of travel times. The slice is along a line through the quiet Sun (shown by the dashed white 
line in the magnetogram of Figure 15) chosen because of the presence of many near-surface large- 
scale flow structures. Inversions for v x at different target depths were performed, the numerical 
divergence d x v x was computed, and the results are given by the colorscale. The color scale is 
such that a positive divergence means an outflow. The v z flows were obtained using f, pi, and 
P2 travel times, and the magnitudes and directions are shown by the arrows. Since we do not 
yet implement a ID depth inversion for v Zl we roughly determine the three depth locations by 
computing the average depth over which the dominant modes of these three ridges probe. We 
emphasize that these placements are only approximate. We also plot the la noise levels of v z at 
each depth. We see that over the whole depth range, the horizontal inflows (outflows) generally 
correspond to vertical downflows (upflows). What one would expect is to see relatively stronger 
vertical flows where the divergence is strongest in absolute value. This is for the most part the 
case. We emphasize that these good correlations are likely not due to cross-talk contamination, 
which tends to diminish as one moves further below the surface (see Section VB). Note also the 
presence of large-scale structures that live for at least 3 days. 
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FIG. 24: Depth slice through a flow map from inversions using 3 days of data. The horizontal divergence of v x 
is shown by the color scale and the arrows represent the velocity of the vertical flows at three different depths. 
The three thick black horizontal lines denote the approximate depth locations of the v z measurements, which 
are found by computing the center of mass of the ID sensitivity functions for the first three ridges (see 
Fig. [9]). The top curve shows v* z , the middle vj 1 , and the bottom v^ 2 . The divergence map was obtained by 
computing maps of v x at several different depths where each depth block, centered about the target depth, 
has a horizontal spatial resolution of < 10 Mm. Then the numerical divergence was calculated. A positive 
divergence means an outflow. The ±la noise for the three v z measurements are given by the two thin lines 
above and below the depth indication line, and all have values less than 10 ms _1 . Note at the top the 
reference scale arrow for v z . 
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IX. SUMMARY AND CONCLUSIONS 

We have presented in detail a fully-consistent procedure for inverting helioseismic travel times to 
infer vector flows in the upper convection zone of the quiet Sun. Travel-time sets are measured for 
all modes that have the same radial order, i.e., along the ridges in the power spectrum. The travel 
times are constructed using an analogue to the common point-to-annulus geometry for 20 annulus 
radii (up to about 30 Mm). Three-dimensional Born sensitivity kernels for the same travel-time 
definition and ridge filtering are computed. In addition, the noise covariance properties of the 
travel times are calculated. Based on the separability of the sensitivity kernels into horizontal and 
vertical components due to the ridge filtering, the inversion is formulated in two steps: the first 
step solves the 2D horizontal problem and the second step solves the ID depth inversion. Optimal 
sets of weights are chosen from both inversions, such that the final averaging kernel is regularized 
in the horizontal and vertical directions. We have provided many examples of averaging kernels, 
which are extremely useful for understanding where in space is the sensitivity of the inversion, as 
well as for determining the amount of cross talk among all of the flow components. It was found 
that the cross talk is reasonably small because the inversion procedure attempts to minimize its 
effect by the use of certain constraints. We furthermore obtain consistent estimates of the noise on 
the measured velocities and the spatial resolution. For practical reasons, the inversion technique 
is convenient since directly inverting for other quantities such as the horizontal divergence of the 
flows or the vertical vorticity only requires one to change the target function. 

Many high-resolution example flow maps have been studied and tested. We have restricted 
ourselves to quiet Sun only. These maps all have horizontal spatial resolution less than about 
11 Mm, or about the pi-mode wavelength at 3 mHz. The recovered flow speeds are below the 
limits for which a linearized theory of travel times is valid [10] . We have tested the inversion in 
several straightforward ways. We have shown that using independent measurements and similar 
averaging kernels gives consistent results. We have also been able to obtain high correlations (~ 0.9) 
with the Doppler velocity data after projecting the inferred horizontal flows onto the line-of-sight 
vector and ignoring pixels with strong magnetic fields. 

We also found that the correlation of 24 h of inferred horizontal flows from day to day on 
average is about 0.4 near the surface and about 0.1 down to about 6 Mm below the surface. 
This is consistent with the conclusion that we are not just measuring noise. However, we find 
that the correlation of flows at a particular depth with the surface flows falls off quite rapidly 
and disappears near ~ 5 Mm beneath the surface, where we do not see any more evidence of 
supergranulation. Similar results have also recently been found on numerical simulations using 
time-distance helioseismology [21] and holography [22]. It could be that for 24 h and at these 
depths the supergranulation signal is completely masked by noise [22l [23] . 

We have shown a direct inversion for the vertical component of the velocity using acoustic and 
surface-gravity waves. The results are in agreement with the overall behavior of the horizontal flows, 
and since the cross talk between v Xl v yi and v z has been made small, we are fairly confident that 
the vertical flows are real. The vertical flows have also been compared to independent inversions 
for the horizontal divergence, and the values are in the expected ranges. We find that the upflow 
speeds in the center of supergranules are approximately 15 — 20% of the horizontal outflow speeds. 
Studying the day-to-day correlations of vertical flow maps also leads us to believe that the signal 
is above the noise. 

Another way to validate many of the findings that we have reported would be to invert the 
available artificial velocity data from realistic numerical simulations of solar convection [26]. Even 
though the averaging kernels give a complete picture of how the data is spatially averaged - a nice 
feature of OLA-type inversions - we intend to carry this out in the near future to study the role 
that noise plays in the interpretation. 
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Of course, we are undertaking many improvements to the inversion presented here. One obvious 
deficiency is the small set of modes we have used. Such a limited number does not allow us to 
obtain many independent target depths, nor any substantially deep ones. More ridges, combined 
with utilizing the spatial frequency content of the waves in each ridge in a more sophisticated way, 
would help us to obtain better, and deeper, averaging kernels. 

Several other improvements currently being studied are ways to minimize the cross talk among 
flow components as much as possible by constructing different types of constraints in the inversion 
procedure. Also, kernels which take into account the line-of-sight projection are almost certain to 
be necessary for inverting data well away from disk center. We already have some of these kernels 
available [27] . 

We thank T. Duvall Jr. for helpful discussions and for providing the data set used in the analysis. 
We also gratefully acknowledge critical comments from a referee that significantly improved this 
paper. SOHO is a collaboration between NASA and ESA. 
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